計算機で一様擬似乱数を生成する．
実数 \(a < b\) に対し，区間 \([a, b]\) 上の 一様分布 (uniform distribution) とは，確率密度関数 \[p(x) = \begin{cases} 0 & \text{if} \ \ x < a \\ \frac{1}{b - a} & \text{if} \ \ a \leq x \leq b \\ 0 & \text{if} \ \ x > b \end{cases}\tag{2.1}\] を持つ確率分布．\(\mathcal{U}[a, b]\) と書く．区間を指定せずに一様分布というと，\(\mathcal{U}[0, 1]\) を指す．
set.seed(1) # seedの設定
runif(3) # 長さ3の一様擬似乱数## [1] 0.2655087 0.3721239 0.5728534関数 runif の出力は一様乱数と認識しても問題ないか？
\(\rightsquigarrow\) スペクトル (Spectral) 検定，コルモゴロフ・スミルノフ (Kolmogorov-Smirnov) 検定 を適用．
# スペクトル検定
set.seed(1)
x <- runif(1000)
y <- runif(1000)
z <- runif(1000)
plot_ly(x=x, y=y, z=z, type="scatter3d", mode="markers", size=0.1)立方体に均等に分布しているように見える．
# コルモゴロフ・スミルノフ検定
set.seed(1)
u <- runif(1000)
ks.test(u, punif)##
## One-sample Kolmogorov-Smirnov test
##
## data: u
## D = 0.024366, p-value = 0.5928
## alternative hypothesis: two-sided一様分布から乖離しているのであれば，\(p\)-値は小さくなるはずであるが，ここでは一様擬似乱数を観測として見たときの \(p\)-値は小さくない．
\(\rightsquigarrow\) 一様擬似乱数 (初期設定はMersenne-Twister法) は一様分布からの独立な乱数列と見ても顕著な差が見られない．
\(a, b, y_0, n\): 事前に決める整数．
\[ y_m = (a y_{m - 1} + b) \mod n \ (m = 1, 2, \ldots) \\ x_m = \frac{y_m}{n}\] として\(\{ x_m: m = 1, 2, \ldots \}\) を出力する．
u <- numeric(3*1e3)
y <- 1234
n <- 2^31 - 1
a <- 65539
b <- 0
for(i in 1:length(u)){
y <- (a*y + b) %% n
u[i] <- y/n
}
u_mat <- matrix(u, nrow=3)
x <- u_mat[1,]
y <- u_mat[2,]
z <- u_mat[3,]
plot_ly(x=x, y=y, z=z, type="scatter3d", mode="markers", size=0.1)2次元平面に整列している様子が見て取れる．
以下，理論的な記述の際は乱数を扱い，実際の計算では擬似乱数を用いるが，擬似乱数も乱数と表記する．
ここで，\(\mathcal{U}[0, 1]\) からの一様乱数が生成できるとする．\(X \sim \mathcal{U}[0, 1]\) に対して，線形変換
\[ Y = b + (a - b) X \] を施すことで，\(\mathcal{U}[a, b]\) に従う乱数が生成できる．
\(\because\) \(\mathcal{U}[a, b]\) は累積分布関数 \[ F(x) = \frac{x - a}{b - a} \ (a \leq x \leq b) \] を持つが，線形変換してできた \(Y\) の従う分布も \[ \mathbb{P}(Y \leq x) = \mathbb{P}\left( X \leq \frac{x - b}{a - b} \right) = \frac{x - b}{a - b} = F(x) \ (a \leq x \leq b) \] となり，同じ累積分布関数を持つため．
一様乱数をもとに，与えられた1次元の確率分布 \(P\) に従う乱数を生成する．
累積分布関数を \(F\) とすると， \[ \mathbb{P}(X \leq x) = F(x) \ (-\infty < x < \infty) \] となるような \(X\) を生成すればよい．
連続な確率分布であれば，任意の \(u \in (0, 1)\) に対して，\(F(x) = u\) となる \(x\) がただ1つに決まる．
\(\rightsquigarrow\) \(F^{-1}(u) = x\) で逆関数 \(F^{-1}\) が定義できる．このとき，一様乱数 \(U\) により， \[ X = F^{-1}(U) \] で \(P\) に従う確率変数を作れる．
\[\because \mathbb{P}(X \leq x) = \mathbb{P}(F^{-1}(U) \leq x) = \mathbb{P}(U \leq F(x)) = F(x) \] このようにして確率分布 \(P\) に従う乱数を生成する方法を 逆変換法 (inversion method) という．
\(\lambda > 0\) とする．指数分布 \(\mathcal{E}(\lambda)\) の累積分布関数は， \[ F(x) = 1 - e^{-\lambda x} \] である．したがって，一様乱数 \(U\) を用いて \[ F^{-1}(U) = - \frac{\log (1 - U)}{\lambda} \] によって指数乱数が発生できる．ここで，\(1 - U \sim \mathcal{U}[0, 1]\) より，\(-\log(U)/\lambda\) によっても指数乱数を生成できる．
R 言語で実装する．パラメータ \(\lambda\) を \(\lambda = 2.0\) とする．
set.seed(1)
lambda <- 2.0
u <- runif(3)
-log(u)/lambda## [1] 0.6630539 0.4942642 0.2785628組み込み関数 rexp でも生成する．
set.seed(1)
lambda <- 2.0
rexp(3, rate=lambda)## [1] 0.37759092 0.59082139 0.07285336同じ seed を用いても異なる値が出力されていることから，異なった実装がなされているようである．
計算時間の比較を行う．
lambda <- 2.0
n <- 1e5
microbenchmark(A <- -log(runif(n))/lambda, times=100)## Unit: milliseconds
## expr min lq mean median uq
## A <- -log(runif(n))/lambda 2.782058 2.837641 2.973486 2.999038 3.043562
## max neval
## 4.524566 100microbenchmark(B <- rexp(n, rate=lambda), times=100)## Unit: milliseconds
## expr min lq mean median uq
## B <- rexp(n, rate = lambda) 3.24646 3.280582 3.397402 3.328264 3.490283
## max neval
## 5.460077 100実行時間に大きな差はないようである．
実数 \(\mu\) と，正の実数 \(\sigma\) に対し，コーシー分布 (Cauchy distribution) \(\mathcal{C}(\mu, \sigma)\) は，確率密度関数 \[ p(x; \mu, \sigma) = \frac{1}{\pi \sigma \left( 1 + \frac{(x - \mu)^2}{\sigma^2} \right)} \ (-\infty < x < \infty) \] を持つ．\(\mu = 0, \sigma = 1\) としたコーシー分布 \(\mathcal{C}(0, 1)\) に注目する．累積分布関数は， \[ F(x) = \int_{-\infty}^x \frac{1}{\pi (1 + y^2)}dy \] となる．学部の解析学を思い出すと，この積分は， \[F(x) = \frac{\mathrm{arctan} x}{\pi} + \frac{1}{2}\] となる．よって， \[F^{-1}(u) = \tan \left( \pi u - \frac{\pi}{2} \right)\] を得る．したがって，一様乱数 \(U\) を用いて， \[X = \tan (\pi U - \pi/2)\] でコーシー乱数を生成できる．
また，コーシー分布は独立な2つの標準正規乱数 \(X_1, X_2\) の比 \[\frac{X_2}{X_1} \tag{2.2}\] の分布と等しい．以下で，逆変換法，(2.2)の方法，組み込み関数の3つの計算時間を比較する．
n <- 1e5
microbenchmark(A <- tan(pi*runif(n)-pi/2), times=100)## Unit: milliseconds
## expr min lq mean median uq
## A <- tan(pi * runif(n) - pi/2) 4.544172 4.667817 4.785406 4.740269 4.863668
## max neval
## 6.416466 100microbenchmark(B <- rnorm(n)/rnorm(n), times=100)## Unit: milliseconds
## expr min lq mean median uq max
## B <- rnorm(n)/rnorm(n) 9.146741 9.350227 9.487853 9.385292 9.489152 12.12001
## neval
## 100microbenchmark(C <- rcauchy(n), times=100)## Unit: milliseconds
## expr min lq mean median uq max neval
## C <- rcauchy(n) 4.731811 4.845782 4.957676 4.916645 5.00278 6.653045 100上述の逆変換法では，累積分布関数 \(F\) の逆関数 \(F^{-1}\) が存在することを仮定した．
離散の場合は逆関数を持たないが，一般に累積分布関数 \(F(x) = \mathbb{P}((-\infty, x])\) は次の性質を持つことを思い出す．
逆に，上記の性質を満たす関数 \(F\) は \(\mathbb{R}\) 上のある確率分布の累積分布関数である．
ここで，関数 \(F\) が単調非減少で右連続のとき， \[F^{-}(u) = \inf \{x \in \mathbb{R}: u \leq F(x) \}\] によって広義の逆関数 \(F^{-}: \mathbb{R} \to \mathbb{R}\) を定める．広義の逆関数 \(F^{-}\) に対しては，\(F^{-}(F(x)) = x\) や \(F(F^{-}(u)) = u\) は必ずしも成り立たないが，次が成り立つ．
\(x \in \mathbb{R}, u \in [0, 1]\) に対して， \[F^{-}(F(x)) \leq x\] \[ F(F^{-}(u)) \geq u \]
確率分布 \(P\) の累積分布関数 \(F\) の広義の逆関数 \(F^{-}\) を用いることで，連続の場合と同様に \(U \sim \mathcal{U}[0, 1]\) から \(F^{-}(U)\) を計算することで，\(P\) に従う乱数を生成できる．
有限集合 \(-\infty < x_1 < \cdots < x_n < \infty\) にのみ値をとる確率分布を考える．この分布の累積分布関数 \(F\) の広義の逆関数 \(F^{-}\) は， \[ F^{-}(u) = \begin{cases} x_1 & \text{if} \ \ 0 \leq u \leq F(x_1) \\ x_2 & \text{if} \ \ F(x_1) < u \leq F(x_2) \\ \vdots \\ x_n & \text{if} \ \ F(x_{n - 1}) < u \leq 1 \end{cases} \] である．可算無限集合でも同様．
幾何分布 \(\mathcal{Ge}(p)\) は，パラメータ \(0 < p < 1\) に対し，確率関数 \[ q(n) = p(1 - p)^n \ (n = 0, 1, \ldots)\] を持つ．累積分布関数は， \[ F(n) = \sum_{m = 0}^n q(m) = \sum_{m = 0}^n p(1 - p)^m = 1 - (1 - p)^{n + 1} \ (n = 0, 1, \ldots) \] である．非負の整数 \(n\) に対し，任意の \(n \leq x < n + 1\) となる実数 \(x\) は \[ F(x) = F(n) = 1 - (1 - p)^{n + 1} \] を満たす．よって，\(x\) の整数部分を \([x]\) と書くと，累積分布関数は \[ F(x) = \begin{cases} 1 - (1 - p)^{[x + 1]} & \text{if} \ x \geq 0 \\ 0 & \text{if} \ x < 0 \end{cases} \] で与えられる．よって，\(F^{-}(u) = x\) となるのは， \[ 1 - (1 - p)^x < u \leq 1 - (1 - p)^{x + 1} \] となるときである．よって，\(F^{-}(u)\) は，\(u > 1 - (1 - p)^x\) を満たす最大の整数なので， \[ F^{-}(u) = [\log (1 - u) / \log (1 - p)] \] と書ける．したがって，\(U \sim \mathcal{U}[0, 1]\) であれば，\([\log (1 - U)/ \log(1 - p)]\) は幾何分布に従う．
p <- 0.2
n <- 3
as.integer(log(runif(n))/log(1 - p)) # as.integer: 実数の整数部分を取り出す## [1] 0 2 6rgeom(n, prob=p)## [1] 4 9 0n <- 1e5
microbenchmark(A <- as.integer(log(runif(n))/log(1 - p)), times=100)## Unit: milliseconds
## expr min lq mean median
## A <- as.integer(log(runif(n))/log(1 - p)) 2.873295 2.896516 2.996405 2.924499
## uq max neval
## 2.992789 4.815373 100microbenchmark(B <- rgeom(n, prob=p), times=100)## Unit: milliseconds
## expr min lq mean median uq max
## B <- rgeom(n, prob = p) 6.684288 6.851178 6.901857 6.894342 6.937812 7.254688
## neval
## 100逆変換法の方が組み込み関数よりも効率的であった．
ポアソン分布 \(\mathcal{P}(\lambda)\) を生成する．幾何分布と同様に，一様乱数 \(U\) および \(n = 0, 1, \ldots\) に対し， \[ F(n - 1) < U \leq F(n) \Longrightarrow X = n \tag{2.3}\] とすれば良い．ただし，形式的に \(F(-1) = 0\) とする．幾何分布と異なり，累積分布関数は有限和 \[ F(n) = \sum_{m = 0}^n \frac{\lambda^n}{m!} e^{-\lambda} \] の形でしか書くことができない．
\(\rightsquigarrow\) 逐次的な条件分岐が必要．\(x = 0, X = x\) として， \[ U > F(x) \Longrightarrow x \leftarrow x + 1 \] なる while ループをする．\(U \leq F(x)\) となったときに while ループを抜け出し，\(X = x\) を出力する．
lambda <- 2
rpoisf <- function(n, lambda) {
c <- exp(-lambda)
res <- numeric(n)
for (i in 1:n) {
x <- 0
q <- c
F <- c
u <- runif(1)
while(u > F) {
x <- x + 1
F <- F + q*lambda/x
}
res[i] <- x
}
return(res)
}
n <- 100
microbenchmark(A <- rpoisf(n, lambda), times=100)## Unit: microseconds
## expr min lq mean median uq max
## A <- rpoisf(n, lambda) 168.919 175.3295 247.7403 178.051 181.8505 6722.139
## neval
## 100microbenchmark(B <- rpois(n, lambda), times=100)## Unit: microseconds
## expr min lq mean median uq max neval
## B <- rpois(n, lambda) 4.065 4.2215 5.24055 4.295 4.505 90.469 100while 文を使うため，組み込み関数 rpois よりもかなり遅い．関数 rpoisf の出力を \(X\) とすると，平均的には \(\mathbb{E}[X + 1] = \lambda + 1\) 回の分岐を行うことになる．
逆変換法の有用性は，逆関数の計算のしやすさで決まる．
\(\rightsquigarrow\) 逆関数の計算コストが高い場合，逆関数を近似で置き換える方法は有効．
ここまでをまとめる．実数上の確率分布 \(P\) の累積分布関数を \(F\) とし，累積分布関数の (一般化された) 逆関数，またはその近似を \(F^{-}\) とすると，以下のようにして \(P\) に従う乱数を生成できる．
変数変換法: 既存の乱数を変数変換して新たな乱数を得る方法
変数変換法は一般の次元に適用できる．以下に変数変換法の考え方を示す．
逆変換法では，一様乱数 \(U\) を関数 \(F^{-1}\) で変換することで，累積分布関数 \(F\) をもつ確率分布に従う乱数を生成した．より一般に，多次元の乱数 \(Y\) を \(S^{-1}\) で変換した \(S^{-1}(Y)\) が目的となる確率分布に従うようにデザインする．
確率変数 \(Y\) が，確率密度関数 \(q(y)\) をもつ確率分布 \(Q\) に従うとする．また，関数 \(S\) が，微分可能で一対一の写像であるとする．さらに，関数 \(S\) の Jacobian を \[ J_S(x) = \text{det} \left[ \frac{\partial S}{\partial x} \right] \] と書く．このとき，確率変数 \(X = S^{-1}(Y)\) に従う確率分布 \(P\) は，正則性の条件の下， \[ p(x) = q(S(x)) |J_S(x)| \tag{A}\] で定義される確率密度関数を持つ．
\(\because\) 変数変換の公式より，開集合 \(\mathcal{O}\) を含む領域で Jacobian matrix が非退化であれば， \[ \mathbb{P}(X \in \mathcal{O}) = \mathbb{P}(S^{-1}(Y) \in \mathcal{O}) \\ = \int_{S^{-1}(y) \in \mathcal{O}} q(y) dy \\ = \int_{x \in \mathcal{O}} q(S(x)) |J_S(x)| dx \] と書くことができる．
よって，(A) が満たされていれば，確率分布 \(P\) に従う乱数を，確率分布 \(Q\) に従う乱数 \(Y\) から生成できる．この事実を利用した乱数生成法が 変数変換法 (general transformation method) である．具体的には，以下のように乱数を生成する．
所望の分布に従う乱数を生成するために，(A) を満たすような \(Q, S\) を見つける必要があるが，見つけ方には自由度がある．
\(\rightsquigarrow\) 以下の性質を利用する．
実数値確率変数 \(Y_1, Y_2\) は独立で，それぞれ確率密度関数 \(q_1, q_2\) をもつ確率分布に従うとする．このとき，\(X_1 = Y_1 + Y_2\) の従う分布は，確率密度関数 \[ q_1 * q_2 (x) = \int q_1(x - y) q_2(y) dy \] を持つ．
\(\because\) 関数 \(S\) を，\(X = (X_1, X_2) = S^{-1}(Y) = (Y_1 + Y_2, Y_2)\) で定める．よって，\(x = (x_1, x_2)\) に対して，\(S(x) = (x_1 - x_2, x_2)\) である．\(S\) の Jacobian は， \[ J_S(x) = \text{det} \left( \begin{array}{cc} \frac{\partial S_1}{\partial x_1} & \frac{\partial S_1}{\partial x_2} \\ \frac{\partial S_2}{\partial x_1} & \frac{\partial S_2}{\partial x_2} \end{array} \right) \\ = \text{det} \left( \begin{array}{cc} 1 & -1 \\ 0 & 1 \end{array} \right) = 1 \] であり，\(Y_1, Y_2\) が独立であることと(A) より，\(X\) の従う分布は， \[ p(x_1, x_2) = q(S(x))|J_S(x)| = q_1(x_1 - x_2) q_2(x_2) \] なる確率密度関数を持つ．したがって，\(x_2\) で積分して，\(X_1 = Y_1 + Y_2\) の従う分布の確率密度関数 \[ p(x_1) = \int_{x_2} p(x_1, x_2) dx_2 = \int q_1(x_1 - x_2) q_2(x_2) dx_2 \] を得る．
以下では，代表的な変数変換法による乱数生成をいくつかまとめる．証明等は省略している．
逆変換法により，指数分布に従う乱数を生成できる．ここでは，指数分布に従う確率変数列を用いて，ガンマ分布，ベータ分布に従う乱数を生成する．
\(\alpha\) は正の実数，\(N\) は正の整数とする．\(X_1, \ldots, X_N \overset{\text{i.i.d.}}{\sim} \mathcal{E}(1)\) とする．このとき， \[ Y = \alpha^{-1} \sum_{n = 1}^N X_n \] はガンマ分布 \(\mathcal{G}(N, \alpha)\) に従う．特に，\(\alpha = 1/2\) とすれば，\(Y\) は自由度 \(2N\) のカイ二乗分布 \(\chi_{2N}^2\) に従う．
よって，逆変換法での指数乱数の生成と合わせると，以下のようにすると \(\mathcal{G}(N, \alpha)\) に従う乱数 \(X\) を生成できる．
alpha <- 2
N <- 10
-1/alpha*log(prod(runif(N)))## [1] 6.894084正の実数 \(p, q, \alpha\) に対し，\(X, Y\) は独立でそれぞれ \(\mathcal{G}(p, \alpha), \mathcal{G}(q, \alpha)\) に従うとする．このとき， \[ Z = \frac{X}{X + Y} \] はベータ分布 \(\mathcal{Be}(p, q)\) に従う．
逆変換法での指数乱数の生成および上のガンマ乱数の生成と合わせると，以下のように \(X \sim \mathcal{Be}(N_1, N_2) \ (N_1, N_2 \in \mathbb{N})\)を生成できる．
n1 <- 10
n2 <- 15
u <- runif(n1 + n2)
sum(log(u[1:n1]))/sum(log(u))## [1] 0.3584183以下の正規乱数生成法を ボックス・ミュラー (Box-Muller) 法 という．
\(U_1, U_2\) は独立で，\(\mathcal{U}[0, 1]\) に従う．このとき， \[ (X_1, X_2) = S^{-1}(U_1, U_2) = (\sqrt{-2 \log U_1} \cos(2 \pi U_2), \sqrt{-2 \log U_1} \sin(2 \pi U_2) ) \] は2次元の標準正規分布に従う．
u <- runif(2)
sqrt(-2*log(u[1]))*c(cos(2*pi*u[2]), sin(2*pi*u[2]))## [1] 1.7109163 0.3909921Box-Muller法のアイデアは極座標表示から来ている．2つの正規乱数 \(X_1, X_2\) をベクトル \(X = (X_1, X_2)\) と捉えると，ベクトルの長さの2乗 \(R^2 = X_1^2 + X_2^2\) は自由度 2 のカイ二乗分布に従う．また，原点から見たベクトルの方向 \(\xi := R^{-1/2} X\) は単位円周上の一様分布に従い，\(R\) と \(\xi\) は独立だと考えられる．
\(\rightsquigarrow\) 逆に，単位円周上の一様分布に従う確率変数として \(\xi = (\cos(2 \pi U_2), \sin (2 \pi U_2))\) を構成し，\(R^2 = -2 \log U_1\) とすれば，\(X = (X_1, X_2) = R^{1/2} \xi\) は独立に標準正規分布に従う．
n <- 1e5
RNGkind(normal.kind = 'default')
microbenchmark(A <- rnorm(n), times=100)## Unit: milliseconds
## expr min lq mean median uq max neval
## A <- rnorm(n) 4.473035 4.592718 4.738723 4.736625 4.841269 6.146236 100RNGkind(normal.kind = 'Box-Muller')
microbenchmark(B <- rnorm(n), times=100)## Unit: milliseconds
## expr min lq mean median uq max neval
## B <- rnorm(n) 3.623579 3.715141 3.841297 3.77388 3.950164 5.482285 100私の環境では，Box-Muller 法の方が逆変換法より計算効率が良かった．
\(X_1, X_2\) は独立な標準正規乱数とする．このとき， \[X_2/X_1\] はコーシー分布 \(\mathcal{C}(0, 1)\) に従う．
X <- rnorm(2)
X[2]/X[1]## [1] 1.936407正の実数 \(\nu\)，標準正規乱数 \(X\), 自由度 \(\nu\) のカイ二乗乱数 \(Y\) に対し， \[ Z = X/\sqrt{Y/\nu} \] は \(\mathcal{T}(0, 1)\) に従う．
nu <- 2.5
X <- rnorm(1)
Y <- rchisq(n=1, df=nu)
X/sqrt(Y/nu)## [1] 0.1333292整数の集合 \(\{1, \ldots, K\}\) に定義された確率分布 \(G\) は確率関数 \(g(k)\) をもつとする．また，\((E, \mathcal{E})\) 上に確率分布 \(F_1, \ldots, F_K\) が定義されているとする．このとき，下記で定まる確率分布を 有限混合分布 (finite mixture distribution) という． \[ P(dx) = \sum_{k = 1}^K g(k) F_k(dx) \]
次の方法で有限混合分布に従う乱数を発生させることができる．
確率分布の列が，正規分布 \(F_k = \mathcal{N}(\mu_k, \sigma_k^2) \ (k = 1, \ldots, K)\) であるとき， \(P\) を 有限混合正規分布 (finite normal mixture distribution) という．確率分布 \(P\) の平均および分散を求める．
\[F_k(dx) = \frac{1}{\sqrt{2 \pi \sigma_k^2}} \exp\left( - \frac{(x - \mu_k)^2}{2} \right) dx \] であるので， \[ P(dx) = \sum_{k = 1}^K g(k) \frac{1}{\sqrt{2 \pi \sigma_k^2}} \exp\left( - \frac{(x - \mu_k)^2}{2} \right) dx \] である．よって，平均は， \[ \bar{\mu} = \int x P(dx) = \int x \sum_{k = 1}^K g(k) \frac{1}{\sqrt{2 \pi \sigma_k^2}} \exp\left( - \frac{(x - \mu_k)^2}{2} \right) dx \\ = \sum_{k = 1}^K \int x\frac{1}{\sqrt{2 \pi \sigma_k^2}} \exp\left( - \frac{(x - \mu_k)^2}{2} \right) dx \\ = \sum_{k = 1}^K g(k) \mu_k \] である．分散は， \[ \int (x - \bar{\mu})^2 P(dx) = \sum_{k = 1}^K \int (x - \bar{\mu})^2 \frac{1}{\sqrt{2 \pi \sigma_k^2}} \exp\left( - \frac{(x - \mu_k)^2}{2} \right) dx \\ = \sum_{k = 1}^K g(k) (\mu_k^2 + \sigma_k^2 - 2 \bar{\mu} \mu_k + \bar{\mu}^2) \\ = \sum_{k = 1}^K g(k) \sigma_k^2 + \sum_{k = 1}^K g(k) (\mu_k - \bar{\mu})^2 \] である．
一般の混合分布を考える．ある状態空間 \((Y, \mathcal{Y})\) 上の確率分布 \(G\) と，\(y \in Y\) をパラメータに持つ \((E, \mathcal{E})\) 上の確率分布の族 \(F_y\) が定義されているとする．また，\(A \in \mathcal{E}\) に対し，\(y \mapsto F_y(A)\) が \(\mathcal{Y}\)-可測とする．
\(\rightsquigarrow\) 混合割合 \(G\) の混合分布 (mixture distribution) が \[ P(dx) = \int_y G(dy) F_y(dx) \] で定義される．混合分布 \(P\) に従う乱数を発生させる手順は，有限混合分布のときと同じである．
パラメータ \((M, \theta) \ (M \in \mathbb{N}, \theta \in (0, 1))\) の負の二項分布 (negative binomial distribution) \(\mathcal{NB}(M, \theta)\) の確率関数は， \[ p(x; M, \theta) = \frac{\Gamma(M + x)}{\Gamma(n)x!} \theta^M (1 - \theta)^x = \binom{M + x - 1}{M - 1} \theta^M (1 - \theta)^x \ (x = 0, 1, \ldots) \] と書くことができる．負の二項分布は， \[ G = \mathcal{G}(M, \theta/(1 - \theta)), \ F_y = \mathcal{P}(y) \] としたときの有限混合分布である．
負の二項分布は，逆変換法をもとにした方法でも生成できる．
確率変数列 \(U_1, \ldots, U_M \overset{i.i.d.}{\sim} \mathcal{U}[0, 1]\) とする．このとき，\(0 \leq \theta \leq 1\) に対し， \[ \sum_{m = 1}^M ([\log(U_m) / \log(1 - \theta)] - 1) \] は負の二項分布 \(\mathcal{NB}(M, \theta)\) に従う．
\(F_Y = \mathcal{N}(0, Y^{-2})\) とする．\(G\) は \(\mathbb{R}_{+}\) 上の確率分布とする．このとき，次で定まる確率分布を，分散混合正規分布 (normal variance mixture distribution) という． \[ P(dx) = \int_{y = 0}^{\infty} \frac{1}{\sqrt{2 \pi y}} \exp \left( - \frac{x^2}{2y} \right) G(dy) dx \]
また，\(P\) に従う乱数は以下の手順で生成される．
ここで，\(t\)-分布は，\(Y \sim \mathcal{G}^{-1}(\nu/2, \nu/2)\) とした分散混合正規分布である．
確率分布 \(P\) が定義されている空間を状態空間 (State space) という．状態空間 \(E\) 上に定義された興味のある分布 \(P\) と同じ状態空間に定義された確率分布 \(Q\) はそれぞれ p.d.f. \(p(x), q(x)\) を持ち， \[r(x) = \frac{p(x)}{q(x)} \leq R \; (x \in E)\] となる \(R > 0\) が存在するとする．
棄却法 (Rejection sampling) は以下の方法で \(Q\) に従う乱数から \(P\) に従う乱数を作り出す．分布 \(Q\) を 提案分布 (Proposal distribution) という．
\(X\) が \(P\) に従うことを示す．簡単のため，1次元の乱数に限定するが，一般の次元で成り立つ．
\(m = 1, 2, \ldots\) について，\(Y_m\) は \(Q\) からの，\(U_m\) は一様分布からの乱数で，すべて独立とする．初めて \(U_m \leq R^{-1} r(Y_m)\) となった \(m \in \mathbb{N}\) を \(\tau\) と書く．このとき，\(X = Y_{\tau}\) とおくと，\(X\) は \(P\) に従う．さらに，\(\mathbb{E}[\tau] = R\).
この定理より，\(P\) に従う乱数を1回生成するのに必要な \(Q\) の乱数の個数 \(\tau\) の期待値は，\(\mathbb{E}[\tau] = R\) である．
\(\rightsquigarrow\) 計算効率的には，期待繰り返し回数 \(R\) が小さければ小さいほど良い．理論上 \[R = \sup_{x \in E} \frac{p(x)}{q(x)}\] とすると最も良い．これは，提案分布 \(Q\) の密度関数が \(P\) の密度関数を「覆う」ことのできる最小の \(R\) となっている．
\(p(\theta)\) を事前分布の事前密度とし，\(p(x | \theta)\) を尤度とする．変数 \(x\) の周辺密度関数 \[p(x) = \int_{\Theta} p(x | \theta) p(\theta) d\theta\] が与えられているとする．\(\hat{\theta}\) を最尤推定量 (Maximum likelihood estimator) とする．
\(\rightsquigarrow\) 事前分布を提案分布として事後分布に従う乱数を生成する棄却法を構成する．
事後密度関数は， \[p(\theta | x) = \frac{p(x | \theta) p(\theta)}{p(x)}\] で与えられるので，密度関数の割合が \[r(\theta) = \frac{p(\theta | x)}{p(\theta)} = \frac{p(x | \theta)}{p(x)}\] となる．\(R\) は関数 \(r(\theta)\) の最大値ととるのが最も効率が良い．最大値は，最尤推定量のときの尤度と周辺密度関数の割合 \[R = \frac{p(x | \hat{\theta})}{p(x)}\] である．よって，この \(R\) を用いて事前分布から事後分布を生成できる．